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Abstract 

We performed molecular dynamics simulations on a one-dimensional diatomic gas 
to investigate the possible long time scale inherent in heterogeneous Hamiltonian 
systems. The exponentially long time scale for energy sharing between the transla- 
tional motion and the vibrational one certainly exists in a large limit of the system 
size. The time scale depends on the vibrational frequency uj not as in the pure 
exponential form ~ exp[£?Lj] but as ~ exp[£?u; a ] with a < 1, in good agreement 
with the expression derived from the Landau- Teller approximation. The numerical 
simulations show that the complete resonance condition for vibrational frequencies 
assumed in the analytical treatment is not essential for the long time scale. Some 
discussions of 1// fluctuations based on the present results will be given. 

Key words: Slow dynamics, MD simulations, Boltzmann-Jeans conjecture, 
thermodynamic limit, heterogeneous systems, l// a fluctuations. 
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1 Introduction 



The experimental fact that the high frequency degrees of freedom do not con- 
tribute to the heat capacity of molecular gases is a seeming violation of the 
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equipartition theorem in the classical statistical mechanics. As a possible ex- 
planation of this effective "freezing" phenomenon, Boltzmann[f] and Jeans 
[2] suggested more than one century ago that, if the time scale in vibration 
and the time scale associated with a typical binary collision in the gas are 
largely different, the amount of the energy exchanged between the vibrational 
and translational degrees of freedom can be small enough so that it takes 
an extremely long time for realization of thermal equilibrium. Although this 
problem was solved by the arrival of quantum mechanics, considerable atten- 
tion has been paid to their original idea (the Boltzmann- Jeans conjecture) as 
a possible scenario of the slow dynamics in heterogeneous systems since their 
viewpoint was first noticed by Benettin et al. [3,4,5]. For example, a modified 
Fermi-Pasta-Ulam model which is composed of two subsystems with differ- 
ent time scales (corresponding to acoustic and optical vibrational modes) was 
analyzed to derive nonequipartition of energy in macroscopic systems [6] . 1/ /- 
type energy fluctuations observed in molecular dynamics (MD) simulations of 
liquid water [7] and the long-term energy storage in protein molecules [8] were 
also discussed in the light of the Boltzmann- Jeans conjecture [9,10]. 

Generic nearly-integrable Hamiltonian systems have the stability time of the 
form 

r N ~ exp [1 /e a ] for < e < e (the Nekhoroshev theorem), (1) 

where e is a positive parameter characterizing the smallness of the perturba- 
tion joined to an unperturbed integrable Hamiltonian; a and e are positive 
constants depending on the explicit form of the Hamiltonian, especially on the 
system size, i.e. the number N of degrees of freedom [11]. The time scale tn 
grows exponentially as the perturbation parameter e becomes small. 

Equation (1) explains the non-ergodic behavior discovered in weakly coupled 
harmonic oscillators by Fermi, Pasta and Ulam (FPU phenomenon) for finite 
N [12]. However, the general proof based on the perturbational techniques 
yields a strong A-dependence for the constants: a ~ 1/N and e ~ N~ N in the 
limit of N — > oo [13], which implies that the FPU phenomenon will disappear 
in large systems or in the thermodynamic limit, although, as pointed out in 
[13], this general results do not exclude the possibility that some specified 
system of the FPU-type could actually exhibit the non-ergodic behavior in 
the large N limit [26]. 

Another type of Hamiltonian systems which can exhibit the long time scale of 
the Nekhoroshev- type is presented in refs. [4,5], the idea being based on the 
Boltzmann- Jeans conjecture. Such systems are represented by 

H = H (p, q) + H^n, + f(p, q, tt, £), (2) 
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where H with the canonical variables p = (pi, • • • ,p u ) and q = (q±, ■ ■ ■ , q v ) 
may generally be non-integrable in contrast to the unperturbed Hamiltonian 
in the Nekhoroshev theorem, while Hi describes the system composed of TV 
harmonic oscillators: 

(3) 




with 7r = (tti, • • • , 71n) and £ = (^i,---, The interaction / is assumed to 
vanish for £ = 0. The stability time is obtained as 

TB~exp[(A/AT], (4) 



where A = uj/VL is the ratio of two typical time scales, I/O. in H and 1/uj in 
H\ and A ^> 1 is assumed. The positive constants a and A* generally depend 
on the number N of harmonic oscillators. 

When the angular frequencies (cui, ■ ■ ■ ,u>n) are nonresonant, the exponent a 
is estimated again as a ~ l/N, but when uji — • ■ • — ujn = (the complete 
resonance) holds, a — 1 is proved for any iV, so that the long time scale of 
the pure exponential form 

t b ~ e A / A * = (5) 



is expected even for macroscopic systems as originally suggested by Jeans. If 
the stability time of the above form remains finite for iV — > oo, the basic idea of 
Boltzmann- Jeans could work as a key to reveal the dynamics of slow relaxation 
observed in macroscopic systems. The situation is, however, not so simple 
because the common perturbational techniques yield strong iV-dependence of 
A* as A* ~ iV 2 even in the complete resonance case [16] and the long time 
scale cannot be guaranteed in the thermodynamic limit. Here too this does 
not mean the long time scale will always disappear for large systems. In fact, 
the above estimate of A^-dependence seems to be statistically unacceptable if 
applied to sufficiently dilute molecular gases in which only binary collisions are 
relevant. Unfortunately, it is not easy to draw definite conclusions about the 
stability time (whether it remains finite or tends to zero in the thermodynamic 
limit) and its more reliable u;-dependence (whether it is pure exponential or 
else) from direct simulations based on the molecular dynamics because the 
computational time steeply increases with the system size and the vibrational 
frequency, though such simulations were attempted for a diatomic gas nearly 
two decades ago [3]. 

Instead of invoking direct simulations, the cu-dependence is considered from 
a different angle in [17]: By assuming a gas composed of identical diatomic 
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molecules (with uj\ = ■ ■ ■ = un = oj) and by applying the Landau- Teller 
approximation [18,19] to binary collisions between molecules the cu-dependence 
of the time scale is obtained as 



exp[Bco a ] with a 



2 2 



(6) 



3 + 2/s 3' 



where the inter-molecular potential 0(r) is of the form 



(j>[r) 



1 



s > 1 for small r. 



(7) 




In the present paper, we perform microcanonical MD simulations of a one- 
dimensional diatonic er dbS SIS db minimal model of macroscopic heterogeneous 
systems (i.e. systems with internal degrees of freedom) and investigate its 
long-term behavior to supplement with some information which has not been 
presented so far, hoping to find a clue to basic understanding of the slow 
relaxation phenomena observed in more realistic systems. 

We first confirm whether the long time scale remains or not in the large N 
limit, which we think is necessary because the time scale does decrease with 
increasing N in a range of small N (for a fixed value of the number den- 
sity). Next, we obtain its cj-dependence to examine the Landau- Teller ap- 
proximation. Since both the perturbation method leading to Eq. (5) and the 
Landau- Teller approximation applied to the binary collisions leading to Eq. 
(6) assume the condition of the complete resonance, we finally test the case 
in which (u>i, • • • , un) are slightly different with each other. 



2 Hamiltonian and Initial Condition 

We consider a one-dimensional molecular gas composed of N identical di- 
atomic molecules with mass M [27]. The positions of two atoms belonging to 
the ith molecule are denoted by ft — a — ^ and ft + a + £i, where ft is the center 
of mass of the molecule and 2a is the bond length, ±£j are the displacements 
of atoms from the equilibrium position. The total Hamiltonian is given by 



H = N [KM + K vih (n) + Kib(0 + Mnt(ft 0} , 



(8) 



where 



1 N v 2 da- 
Ktr (p) = -Y— with Pl = M^- 



(9) 
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is the kinetic energy of translational motion of molecules, 

1 N 7T 2 df- 
K vib (7l) = — V-^r With 7Tj = M~y~ (10) 



and 



^ ( «4g^ (id 



are the kinetic and potential energies of intra-molecular vibrations, respec- 
tively, and 

VW<Z, = 4 f> ((ft - a - CO " (ft-i + « + 6-0) (12) 
iV i=i 



is the inter-molecular potential. We employ the cyclic boundary condition: 
g = q N and £o = 6v- The molecular interaction is identical with that adopted 
in [3], which is short range repulsive and given by 

e -(r/ro) 2 

0(r) = , r > 0. (13) 

r/r 



In the present system, N [K tT (p) + V int (q, 0)] corresponds to H (p, q) with v = 
N and JV[Mnt(g,0 -MntM)] to /(p,g,7r,0 in Eq. (2). We may set M = 
r o = 0o = 1 without losing generality by choosing the units of mass, length 
and energy as M, tq and <f>o, respectively. Then the unit of time is ro^M/00 = 
2tt/Q, which is a typical time scale characterizing the Hamiltonian Ho(p,q). 
Note that u scaled by this time scale is actually read as the time scale ratio 
u/Q of the two subsystems. 

Not like many FPU-type simulations starting from a state in which the energy 
is distributed to only specific normal modes, we chose a nearly thermal equi- 
librium state as the initial state, except in Section 3.5. This is because we are 
mainly interested in the long-term behavior observed in macroscopic systems 
in thermal equilibrium. 

If the system represented by the Hamiltonian (8) is in a thermal equilibrium 
state with the temperature T, the expectation value of the kinetic energy K tI 
(in the canonical ensemble) is calculated as 



<*tr(p)> = 



j? S£i J Ae- H / kT dqdtdpdn 
J e- H / kT dqd£dpdTr 
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If /^e-rf/^cfo kT 

N^i f e -tf/ 2MkT d Pl 2 • 1 ; 

Note that Eq. (14) is exact even if the interaction term Vj n t((?, £) exists, because 
the variable is separated from the rest of the variables. Similarly, one obtains 

kT 

(KM) = \, (15) 
which is also exact. 

Without the term Vj nt , the Hamiltonian is decoulped into three terms with 
N degrees of freedom for each: NK tT (p), NK v fo(ir), and iWvib(0> an d the 
expectation value of Kdb(0 i n the thermal equilibrium state is also obtained 
rigorously as 

kT 

(Kib(O) = T , (16) 



namely , the energy is equally distributed among three "subsystems" as 

(NK tT (p)) = (NK vih (n)) = (iVKib(O) = (17) 



The present system is isolated (not contacted with the heat bath) and includes 
the coupling term Vint- However, if the coupling energy is relatively small and 
iV is large, the energy of each subsystem in the stationary state is expected 
to be approximately equal, namely the long-time average for each subsystem 
satisfies 



NEi 



NK tI (p) ~ NK vih (n) ~ W vib (0 ~ (18) 



where E tot is the specific energy, i.e. the total energy of the system per 
molecule. 

Taking account of the above, we start simulations with an initial state which 
is approximately thermal equilibrium as follows: 

For the coordinate £, we generate iV independent random numbers R^^ (i = 
1,2, ■ ■ ■ , N) with the normal distribution (Gaussian distribution with zero 
mean and unit variance) and set the initial value as £j(0) = S^R^i, where 
the positive constant is chosen so that 

_ i f M^goy _ Etot 

Kib ' 0= iV^ 2 -"3" (19) 
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holds. This condition corresponds to Eq. (18), i.e. a thermal equilibrium state. 



The initial position of the center of mass for each molecule is chosen as 

<&(()) = (i-l)L + 5,^, (20) 



where L is the average molecular distance, S q R q ^ is a small deviation from 
the average and R q ^ (i — 1, • • • , N) are N independent random numbers with 
the normal distribution. The constant S q is a small positive number so that 
the condition ^(0) < q i+ i(0) (i — 1, • • • , N) is satisfied. The initial value of 
the molecular interaction energy 

vW> = Mnt(g(0U(0)) (2i) 



is expected to be small by choosing the inter-molecular distance L being large 
enough in comparison to the potential parameter r . 

Next, we choose the initial momentum Pi(0) (i = 1,---,N) by a similar 
way as for £j(0) but with a constraint X^LiPi(O) = 0: We generate N in- 
dependent random numbers R p ^ with the normal distribution [28] and set 
Pi(0) = S p (R p ,i — J2f=i Rp,j/N^j with a constant S p which satisfies 

^• = iV^W = 2 • (22) 

i=i 



Then K tr is almost equal to E tot /3 when V^o is small. 

Finally 7Tj(0) (i = 1, 2, • • • , N) are chosen from N independent random num- 
bers R^ ti with normal distribution and set 7Tj(0) = S w R n:i with a constant Un- 
satisfying 

K vih , = J^J2 = K tr,o- (23) 

i=l 



Hereafter we mean the initial state satisfying Eqs. (19) ~ (23) by "nearly 
equilibrium initial state". Since the parameters L and a appear only in the 
form L — 2a in Vint ( see Eq. (12)), we may set a = by including its effect in 
the parameter L. 
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Fig. 1. A sample of the temporal fluctuations of energies for iV = 128 
and uj = 8. E tot : the total energy (per molecule), i^trF the ki- 
netic energy of the translational motion, = K v[ ^ + V V ^F the vibra- 
tional energy, VJntFthe molecular interaction energyD The initial condition: 
Etot,o = l,Kib,o = l/3,T4nt,o = 0.000 024,K trj0 = ^ vibj0 = 0.333 322. Instanta- 
neous values are plotted with the sampling time interval At = 0.04 in the interval 
[230, 230 + 40] after the system was relaxed for i rc i ax = 230. 

3 Numerical Results 



The computer simulations were performed using the symplectic method [20]. 
As mentioned in the previous section, We set M — ro — 0o — 1- Throughout 
the present simulations the average molecular distance was fixed as L = 4. 



3.1 Time evolution 



Employing the initial condition in Eqs. (19) ~ (23), we started the simula- 
tion from a nearly thermal equilibrium state, in which K tl (p) = K v i^(7i) m 
Knb(0 — Etot/ 3 an d Vi n t ~ 0. The time dependence of energies of the total 
system (per molecule) for a short time interval 230 <t< 230 + 40 is presented 
in Fig. 1 for AT = 128 and uj — 8, in which 

E tot = K tT + E vih + (24) 



and 



E. 



vib 
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vib 
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vib' 



(25) 



8 



1.0 



E, 



tot 



0.8 ■ 



CD 
LU 




Fig. 2. Long time behavior. The same system parameters and the same initial 
condition as in Fig. 1. Instantaneous values are plotted at 4096 time points with 
the sampling interval At = 1500. 

Figures 2, 3 and 4 show the long time behavior of energies and their standard 
deviations. In Fig. 3, the plotted point at the time t is the average of the 
quantity A defined by 



In Fig. 4, the standard deviation at the time t for the quantity A is plotted, 
which is defined by 



Each energy and its standard deviation are plotted with the sampling time in- 
terval At = 1500. The final values of energies in Fig. 3 are E tot = 1.000 001 63, 
K tr = 0.310 09, E vih = 0.312 39 + 0.311 86 = 0.624 25, V int = 0.065 66. Those 
values agree well with the thermal equilibrium energies K tr = (E tot — Vi nt )(N — 
l)/(3iV - 1) = 0.309 82, K vih = V vih = (E tot - V int )N/(3N - 1) = 0.312 26 
for N = 128. The factors (N - l)/(3iV - 1) and N/(3N - 1) instead of 1/3 
appear because the center of mass of the total system is at rest. In Fig. 4, the 
final value of the standard deviation a of the total energy E tot is 1.2 x 10~ 5 , 
which is a measure of the computational error. 

It should be noted that the energies of individual molecules change quite 
largely (compare with E tot = 1 which is averagely distributed to each molecule) 
as can be seen in Fig. 5 in which instantaneous energies of molecule 1 are plot- 
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Fig. 3. Long time behavior of energy fluctuations. The same run as in Fig. 2, but 
the plotted point at the time t represents the average over the time interval [0, t] 
defined by (26). The sampling interval is the same as in Fig. 2. 



0.05 




Fig. 4. Standard deviation a of energies corresponding to Fig. 3. Plotted points at 
the time t are obtained from (27). The sampling interval is the same as in Fig. 2. 
The standard deviation of the total energy, a of E tot , is 1.2 x 10 -5 . 



ted. The large energy exchange is induced by the collision with the neighboring 
molecules as seen from Fig. 6. 
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Fig. 5. Energies of molecule 1 in the same run as in Fig. 1. Ki jtr = p\/2M, 
E iyih = irf/2M + Moj 2 g/2, V 1M = {4>{qi -q -£i- &) + Hl2 ~ qi '- 6 - £i))A 
^l.tot = Kl,tr + ^l.vib + ^l,int- 
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Fig. 6. Motion of atoms belonging to molecules 0, 1 and 2. The coordinates of 
atoms qi ± £j (i = 0, 1, 2) are plotted in the same run as in Fig. 1. 

3.2 Size- dependence of the correlation time 



As a time scale characterizing the system in equilibrium, one may employ 
the correlation time of the energy fluctuation. Following [3] , we computed the 
correlation time r cor of the vibrational energy E vih which is defined such that 
the normalized correlation function 



G(t) = 



E vih (t)E vih (t + r) - (E vih (t)y 
E^^-(E~^tj) 2 



(28) 
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Fig. 7. Correlation function G(t) of the vibrational energy E v ^ for 20 independent 
runs with nearly equilibrium initial states. N = 128, to = 8, £tot = 1- 

becomes a half of the initial value G(0) = 1: 



G(r cor ) = 0.5, 



(29) 



where the overline denotes the time average over the dummy variable t for a 
sufficiently large interval t Q b s - 

We computed r cor for various samples with different random initial conditions 
satisfying Eqs. (19) ~ (23) and define the average by 



Ave [r cor ] . 



(30) 



We have also employed another averaged correlation time t' cot which is defined 
such that 



G(0 = 0.5, 



(31) 



where G(r) is the correlation function averaged over the samples: 

G(r) = Ave [G(r)} . (32) 



The correlation time t' cox is more easily estimated especially when G(r) decays 
very slowly in some samples and it takes very long time to obtain the average 
r cor over all samples, which actually occurs for small N. 

The correlation function G(r) is plotted in Fig. 7 for iV = 128, cu — 8, E tot = 1 
and for 20 samples with different random initial states satisfying (19) ~ (23). 
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Fig. 8. System size dependence of the t ime scale. Two kinds of correlation time T CO r 
defined by (30) and t' cov by eq. (31) are plotted against the number N of molecules. 
The average is taken over 20 samples, i.e. 20 independent runs starting from nearly 
equilibrium initial states. The error bars indicate r cor ± a, where a is the standard 
deviation for 20 samples, u = 8, E tot = 1.0. 

To estimate G(t) reliably over the range < r < 300, the time average was 
taken in the interval [t re iax, ^rciax + ^obs] with the observation time t G bs = 1 x 10 6 
after the system was relaxed for t rclax = 2 x 10 4 . 

In Fig. 8, r cor and r' COT are plotted against the system size, i.e. the number N 
of molecules, from N = 4 to N = 8192, where tv = 8 and the specific energy is 
fixed at E tot = 1.0. The average is obtained from 20 independent runs starting 
from nearly equilibrium initial states and the error bars indicate ±er, a being 
the standard deviation of r cor for 20 runs.. The correlation time as a function 
of N decreases with increasing N in the beginning but tends to a finite value 
pa 71 as iV — * oo. The difference between r cor and r' cor is large for small values 
of N but very small for large N. The error bar of r cor also becomes small for 
large iV if the precision of the numerical simulation is high enough. Figure 8 
proves that the correlation time certainly stays finite for large N, where the 
molecular density is fixed at l/L = 1/4. 

We have obtained a similar size dependence for a larger value of the specific 
energy: When E tot = 2.0 with u = 8, both r cor and r^ or tend to ~ 10.9 with 
increasing iV as (N,T cor ± a,r^ or ) = (8, 13.6 ± 4.7, 13.0), (16, 11.0 ± 1.0, 10.9), 
(32, 11.1 ± 0.6, 11.0), (64, 11.0 ± 0.3, 11.0), (128, 10.9 ± 0.2, 11.0), (256, 10.9 ± 
0.1, 11.0), (512, 10.9 ± 0.1, 10.9). For larger values of E tot , naturally the time 
scale is smaller but reaches a constant value faster. Similar dependence of the 
time scale on the system size and on the specific energy has been reported 
from the numerical simulations of FPU /3-model [14]. 
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03 Fig.9 

Fig. 9. Semilog plot of w-dependence of the averaged correlation times T cor and r' cor . 
The average is taken over 20 independent runs. Error bars indicate T cor ± cr, where 
a is the standard deviation for 20 runs with nearly equilibrium initial states. The 
fitting lines are r cor = 0.0090 x exp[3.92 x uj a ) and t' cot = 0.0096 x exp[3.88 x a; - 4 ] 
with a = 0.4. N = 128, E tot = 1.0. 

3.3 u- dependence of the correlation time 



In Fig. 9, the averaged correlation times r cor and t' cot are plotted against u for 
N = 128. If the mathematical theorem (5) is applied directly to the present 
system in which the condition oj\ = ■ ■ ■ = ojn = 0J is satisfied, the plot should 
drop on the straight line: ~ exp[Bu a ], with a — 1, while the present results 
can be fitted better with a smaller than one. The lines are for a = 0.4, which 
was suggested from Eq. (6) (the Landau- Teller approximation) with s = 1 
corresponding to lim r ^ o 0( r ) ~ ^/ r ' s f° r the potential </>(r) in Eq. (13). We 
obtained the fitting curves r cor = A exp[Bcu 0A ] with (A, B) = (0.0090,3.92) 
and t' cot = A' exp{B'tu 0A } with (A', B ) = (0.0096,3.88) [30]. Also, Fig. 10 
shows good fitting to the same ^-dependence in case of iV = 256 with (A, B) = 
(0.0096, 3.88) for r cor and (A', B') = (0.0102, 3.85) for t' cot . 

Deviation of the (^-dependence from the pure exponential form can be qual- 
itatively understood as follows [17]: The energy transfer AE from the trans- 
lational motion to the vibrational one arises dominantly from direct binary 
collisions of molecules. The two subsystems of translational and vibrational 
degrees of freedom are assumed to be in equilibrium states with well-defined 
temperatures T tr and T vib , respectively during the single binary collision. The 
Landau- Teller approximation applied to a single binary collision yields AE as 
a function of the asymptotic data of two molecules before the collision, which 
is of the form AE ~ e~ TUJ . The parameter r depends on the asymptotic value 
p (at t — — oo) of the relative momentum of the molecules, or equivalently 
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Fig. 10. The same as in Fig. 9 but for N = 256. The fitting lines are 
r cor = 0.0096 x exp[3.88 x uj a ] and t' cov = 0.102 x exp[3.85 x u a ] with a = 0.4. 



on the translational energy E tr ~ p 2 /2, and is expressed as r ~ (_E tr )~ 7 with 
7 = (s + 2)/2s > for the inter-molecular potential (7). The average energy 
exchange per unit time is given by integrating AE over various asymptotic 
states with the Boltzmann factor exp[— E tT /kT tr — E vih /kT vih \. Then the final 
form of the main w-dependence of the average energy exchange per unit time 
becomes stretched exponential: 



(AE) ~ J exp [-(£ tr )-^] exp [-£ tr /A;T tr ] d£ tr 



with 



(33) 



■exp[-Bcu a ] 



a = 1 — 



7 



1 + 7 3 + 2/s' 



(34) 



which leads to the time scale ~ 1/ (AE) as in Eq. (6). The origin of the expo- 
nent a smaller than one arises from the balance of the negative exponent —7 
in the parameter r, (i.e. the energy exchange AE ~ e - ™ in individual colli- 
sions increasing with E tT ) and the weight of the Boltzmann factor decreasing 
with E tT . 
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Iog 10 f log 10 f 

Fig. 11. (a) Log- log plot of PSD for the vibrational energy: and for translational 
energy: K tr . The spectra for energies of a single molecule are denoted by vib an d 
Ki,tv N = 128, oj = 16, -Etot = 1-0. The time interval of FFT time series data is 
At = 0.05, the number of the FFT data is iVdat = 2 26 , i.e. the observation time is 
t bs = 3.4 x 10 6 . Each PSD was obtained from the average over 20 independent runs 
with nearly equilibrium initial states, (b) The same as in (a) but with the initial 
states which are the final states of 20 runs in (a). 

3.4 Power spectra and various time scales 

The power spectrum density (PSD) S(f) of energy fluctuations is presented 
in Fig. 11(a) for vibrational motion and for translational one, where oj = 16, 
N = 128, E tot = 1.0 and the observation time t Q b s = 3.4 x 10 6 . The sharp 
peak at / v ib = oj/2tt m 2.5 in the spectrum of E vi ^ corresponds to the smallest 
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Fig. 12. Semilog plot of the distribution function n{x) of molecular distance 
for N = 128,1024 and N = 8192. u = 8,E tot = 1.0. n(x) is normalized as 
f n(x) dx = N . The initial distribution [t = 0) and the distribution at t = 200 
are compared. At t = 0, the distribution is gaussian with the center on x = 4. The 
distribution at t = 200 is almost final one and is well fitted to ~ e ~ x / x o with x ~ 3 
for x ^> xq. Each line was obtained from the average over 200 independent runs 
with nearly equilibrium initial states. 

time scale of the present system: t vib = l// v ib ~ 0.39. (The second peak is its 
harmonics at / = 2/ vi b.) 

The next small time scale (apart from the time scale of the inter-molecular 
potential: t un i t = r ^jM/(pQ = 1) is the mean free time (collision time): t co i = 
L/v ~ 4.9, where the molecular distance is L = 4 and the average velocity of 
the molecule is v 0.82 because the average translational energy is Mv 2 /2 = 
K tT ps 1/3. Corresponding to this time scale, a mild hump around / co i = 
1/^coi ~ 0.21 can be seen in the PSD of K tr . 

We have employed the initial states in which the positions (i — 1, 2, • • • N) 
of molecules are distributed with Gaussian distribution with a small width 
around the equilibrium position as in (20). The distribution function n(x) 
of the molecular distance x = pj+i — qi is expected to tend to exponential: 
n(x) ~ e ~ x / x o f or large x after a relatively short time which is long enough 
in comparison to t col . This is because the molecular interaction is short range 
repulsive and each molecule occupies its position on a line independently. 
That is confirmed by semilog plot of n(x) for N = 128, 1024 and N = 8182 
in Figure 12. The initial gaussian distribution with the center x = L tends to 
exponential with the average value Xq ~ 3 as can be explained in terms of the 
average distance L = 4 and the scale of the repulsive potential r — 1 which 
yields Xq ~ L — r = 3. 
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Fig. 13. The same as Fig. 11 but with a shorter observation time: t bs — 

1.6 x 10 4 . 

At = 0.5 and A^at = 2 15 . The spectra of and -fT v ib is n °t yet saturated and 
look power-low like, while the spectra for the single molecule (-Ei iV ib an d -?^i,vib) ar e 
nearly flat in a low frequency range. 

The spectra of E v ^ and K tr are flat at frequencies lower than a certain fre- 
quency at the "shoulder" which we denote f . In Fig. 11(a), the shoulder 
frequency is / ~ 10" 5 . To observe the saturated spectrum (the spectrum 
which is white at low frequencies) one needs to observe the system for a time 
interval at least to = l//o ~ 10 5 . In other words, to is a typical time scale after 
which the system loses the memory of the initial state. It should be pointed 
out that this long time scale characterizes the system in thermal equilibrium 
and does not depend on the specified initial state if chosen as an equilibrium 
one. To confirm it, we performed the simulation starting from the final state at 
t = tobs of each run in (a) of Fig. 11 and continued the simulation for another 
time interval t G b s - The power spectra thus obtained were shown in (b) of Fig. 
11, which is very similar to those in (a). 

We have defined r cor as a characteristic time scale of the system, which is 
r cor pa 1.1 x 10 3 in the case of N = 128, u — 16 (see Fig. 9). If the correlation 
function is pure exponential (Debye-type relaxation) as G(t) ~ e~ 27T - fcT , the 
corresponding PSD is Lorentzian given by S(f) ~ 1/[1 + (f / f c ) 2 ] from the 
Wiener-Khinchin theorem and S(f) is flat at low frequencies / < f c , i.e. 
f c is the shoulder frequency. If our G(t) were exponential, this frequency 
should relate to the correlation time T cor defined by (29) as f c — ln2/27TT co i ~ 
ln2/27rr co i ps 1.0 x 1CT 4 . As we have seen in the above, the shoulder frequency 
/o is much smaller (at least more than one digit) than f c , which means that 
decay of the correlation function is slower than exponential, i.e. the PSD 
function deviates from Lorentzian. (See below.) 

The long time scale appears clearly when one observes the energy involved with 
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all molecules but not the energy of individual ones. The spectra of the vibra- 
tional and translational energies of molecule 1 (£i )V ib = ^i/2M + Mu 2 ^\/2 
and K 1>tT = pf/2M) are also presented in Fig. 11. The shoulder frequencies 
in these spectra are much higher (order of three or four digits) than in the 
spectra of E v ^ and K vi t,, implying that each molecule loses its initial memory 
much faster than the total system. Similar results for the energy power spectra 
of the total system and the individual molecules have been obtained from MD 
simulations of liquid water [9]. 

The long time scale is often related to the power law of PSD as ~ 1/ f a with 
< a < 2 over a wide range of frequency. It should be noted, however, that 
the spectrum sometimes appears to be 1 / /"-like when the observation time 
tobs is not long enough. Power spectra obtained from the same simulations as 
in Fig 11 but with a shorter observation time: t obs ~ 1.6 x 10 4 and a longer 
sampling time At = 0.5 are presented in Fig. 13 [29]. For energies of the 
single molecule, the power spectra are saturated, i.e. white at low frequencies, 
while the spectra of E v ^ and K tr look like ~ l// a with a being 1.4 ~ 1.5 
over a low frequency range of about two digits and no shoulders are observed 
because the lowest frequency l/t Q bs is higher than f . As mentioned above, 
these spectra deviate from the Lorentzian form in which a = 2 should be 
expected for / > /q. The value of a smaller than two might be attributed 
to the heterogeneity, i.e. the system is more complex because of the internal 
degrees of freedom. In fact, MD simulations of much more complex systems 
[9] (three-dimensional liquid water, alcohol and argon cluster) yielded power 
spectra whose exponent a of the power-law decaying part is smaller than the 
present value. 

3. 5 Relaxation from nonequilibrium initial states 

In deriving the time scale in Eq. (6) from the Landau- Teller approximation, a 
nonequilibrium initial state is assumed and the relaxation time toward equi- 
librium is estimated. Here we perform simulations based on this idea and 
initially prepare the system with different temperatures for subsystems of the 
translational and vibrational motion. We set the energies of subsystems as 
Ktr,o = 0.2 and -E v ib,o = -^vib,o + Kib,o — 0.8 corresponding to the initial 
temperatures kT trfi = 2K tr ^ = 0.4 and A;T vibi0 = -E v ib,o = 0.8. In Fig. 14, the 
temporal variation of the quantity E vih — 2K tr = k(T vih — T tr ) is presented 
for various values of u, which shows how the temperature difference of two 
subsystems tends to zero and the thermal equilibrium is realized. 

If the temperature difference AT = T vib — T tr decays as 

d .. AT(t) , . 

—AT(t) = (35) 
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Fig. 14. Temperature difference between the vibrational and translational mo- 
tion. N = 128. Temporal variation of kAT = kT vih - kT tr = E vih - 2K tI is 
plotted for to = 8, 12, 16. Each line was obtained from the average over 200 
independent runs with different initial states satisfying the condition Etot.o = 1, 
Kvib,o = Kib,o = 0.4, K tI fi = 0.2, Vi nt ,o ~ so that the initial temperature is 
&^vib,o = 0.8, kT tTt o = 0.4. The plotted point at the time t represents the average 
over the time interval [t — 100, t+ 100] . The sampling interval of plotting is At = 100. 
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Fig. 15. Semilog plot of w-dependence of the relaxation time t\ defined in Eq. (35). 
N = 128, Etot = 1.0. The same initial condition as in Fig. 14. The guide line is 
A 1 exp [Siw - 4 ] with A 1 = 0.204 and B x = 3.601. 

the relaxation time T\ derived from the Landau- Teller approximation depends 
on uj and T tr , and its main ^-dependence is obtained as T\ = C exp[Buj a ], 
where C still depends on uj weakly and on T tr [17]. If T tr -dependence of C could 
be ignored, the lines in Fig. 14 should be fitted to fcAT(i) = £;AT(0)e-*/ T1 
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with kAT(0) = 0.4. In fact, these lines exhibit some oscillations as a function 
of the time t, even though each line has been obtained from the average over 
the time interval [100 — t,t+ 100] and from the average over 200 independent 
runs. Especially the lines (for larger values of u) starting from 0.4 do not 
decrease in the beginning but increase steeply (up to ps 0.47 for uo = 16) [31]. 
Hence, the fitting to the exponentially decreasing function is very rough. We, 
however, estimated the relaxation time T\ using the fitting function 0.4 e - */ 71 . 
The a>-dependence of T\ is given in Fig. 15 and the guide line is obtained from 
T\ = A 1 exp [i^u; ' 4 ] with A 1 = 0.240 and B x = 3.601. In comparison to the 
correlation time (see Fig. 9), data points of T\ appreciably deviate from the 
fitting line. The main reason might be the temperature dependence of t\ in 
Eq. (35) which has been ignored. The correlation time is estimated in thermal 
equilibrium, while the relaxation time here is defined in the nonequilibrium 
state, and the smaller the temperature difference, the better the approximation 
[32]. 



Next we compare the relaxation time T\ (the time scale of energy exchange 
between the vibrational and the translational motion) and another time scale: 
the time scale of energy relaxation within the vibrational degrees of freedom, 
we carried out a simulation starting from a nonequilibrium vibrational state, 
i.e. the initial state is such that the vibrational motion of subsystems I (com- 
posed of the first N/2 molecules) and II (composed of the rest N/2 molecules) 
are in different thermal equilibrium, while the translational motion is in an- 
other thermal equilibrium state as a whole. The initial values of vibrational 
energies (per molecule) -Ei, v ib, -En.vib of subsystems I and II were chosen as 
#i,vib,o = 2£ vibi0 and £n, v ib,o = 0, respectively. 



Fig. 16 shows the case for N = 256, u — 16 and for the initial condition 
E v ib,o — 0.8 and i£tr,o — 0.2, where plotted points were obtained by the average 
over 20 different runs. The temperature difference defined by kAT = E v ^ — 
2K tr , starting from kAT = 0.4, decays monotonically, while the temperature 
difference kAT v ^ = — -En )V ib exhibits a damping oscillation. For a fixed 
value of uj, the first zero-crossing time of kAT v ^ increases linearly with N, and 
the amplitude of the damping oscillation, i.e. the most negative value of kAT vl ^ 
turned out to decrease with increasing N. On the other hand, for a fixed value 
of N, the first zero-crossing time of kAT vih increases weakly with increasing ui 
but Ti, i.e. the time scale of kAT increases much strongly (~ exp Bw ° A ). Thus, 
for sufficiently large values of both N and u, the time scale of energy exchange 
between the translational and vibrational motion is much larger than that of 
energy relaxation within the vibrational system, which can actually be seen 
in Fig. 16. 
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Fig. 16. Temporal variation of the quantities kAT = E v n, — 2K tr and 
/cAT v ib = £^i, v ib ~~ ^n,vib) where the suffixes I and II indicate the subsystems com- 
posed of the first N/2 and the rest N/2 molecules, respectively. Initial condition: 
-Etot,o = 1) -Evib,o = 0.8, K tr fi = 0.2, Vjnt.o ~ 0. The vibrational energy is distributed 
only to the subsystem I so that i^i )V ib,o = ^i,vib,o = 0.8, and i^n iV ib,o = ^n,vib,o = 0. 
Each line was obtained from the average over 20 independent runs satisfying the 
same initial condition. N = 256, id = 16. The sampling interval of plotting is At = 5. 

3. 6 Effect of deviation from the complete resonance 

As mentioned already, the exponent a in Eq. (4) generally depends on the num- 
ber N of degrees of freedom, and the long time scale might not be realized for 
large N except when all u^s are equal to each other. Also in the Landau- Teller 
approximation, the binary collisions are such that the vibrational frequency 
of each molecule is identical. Hence it is interesting to investigate whether the 
condition of the complete resonance is essential for the long time scale in the 
system with many degrees of freedom. 

We have chosen homogeneous random numbers (cui, • • • , u>n) in the interval 
[uj — Auj/2,uj + Aid/ 2], whose average is id. Fig. 17 shows Atu-dependence 
of the correlation time for N = 128, E tot — 1, u — 8. The correlation time 
decreased with increasing Au; monotonically, but its Au;-dependence is rather 
complicated, which we have not yet investigated. 

By assuming a fixed small value for Aid as Aid = 0.1, we obtained the in- 
dependence of the correlation time as in Fig. 18. The present result shows 
that the condition of complete resonance: id\ — ■ • • — is not essential for 
the long time scale because the correlation time increases with id in the same 
way as in Aid = for a finite (but small) value of Aid. 
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Fig. 17. Aw-dependence of the averaged correlation time r cor . N = 128, E tot = 1.0. 
The vibrational frequencies of molecules ■ ■ ■ , wjv) are generated from homoge- 
neous random numbers with the width Aw and a fixed average value lo = 8. r cor is 
obtained from the average over 40 independent runs with nearly equilibrium initial 
states. Error bars indicate r cor ± cr, where a is the standard deviation. The points 
of Alj = correspond to those of N = 128 in Fig. 9 but obtained from 40 runs. 




Fig. 18. Semilog plot of w-dependence of r cor for a fixed value of Aw = 0.1. 
N = 128, -Etot = 1-0. The average is taken over 40 independent runs with nearly 
equilibrium initial states and error bars express ±a. The guide line corresponds to 
r CO r = Aexp[B x lo oa ] with A = 0.0106, B = 3.83. 

4 Conclusions and Discussions 

In conclusion, we estimated the correlation time r cor of the vibrational energy 
which is denned in thermal equilibrium. The present simulations yielded the 
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time scale which indeed remains finite in the large N limit and grows as 
~ expi?u; a4 in good agreement with that obtained from the Landau- Teller 
approximation. The condition of the complete resonance assumed both in the 
mathematical theorem based on the perturbation method and in the Landau- 
Teller approximation is not essential for the long time scale because a slight 
deviation from the complete resonance did not change the tu-dependence of 
the time scale. 

Concerning the system size dependence, we found that the correlation time 
T cor(N) decreased in the beginning but tends to a finite value with increasing 
N for a fixed value of the molecular density. If the molecular density p = 
N/Ltot = 1/L (L to t is the length of the system) is small enough so that 
three-molecule collisions can be neglected, then one can expect r cor (iV) being 
constant (for a fixed value of p) . In other words, the system size dependence 
of the time scale is regarded as arising from multi-molecular collisions. Then 
a natural question arises: Does the long time scale in the thermodynamic 
limit remains for high density case? The present paper does not answer this 
question. 

The present study was done for a one-dimensional gas. To what extent does 
the one-dimensional scenario work in higher dimensional gases? On this issue, 
we note the case of FPU-type systems (as a model of lattice vibrations) [26]: 
The long time scale remains in a one-dimensional chain but disappeared in a 
two-dimensional lattice, which suggests the difficulty in predicting a correct 
answer for the gas model too. No systematic study of the size dependence 
has been reported for higher-dimensional systems including solids and liquid. 
Nevertheless, our feeling is that, as suggested in [9], if the system is heteroge- 
neous (i.e. composed of subsystems with well separated time scales), the basic 
idea of Boltzmann and Jeans might be still usable and so the slow relaxation 
might be realizable in the thermodynamic limit for the high-dimensional gases 
as well and even for more general systems covering liquid and solid. 

The power spectra of the energy fluctuations obtained in the present simu- 
lations are 1 / /"-like with a around 1.4 ~ 1.5 in a range of low frequencies. 
In "real" 1/ / fluctuations (or noise) observed in, say, the electrical resistance 
of condensed matter systems, a is nearly constant and very close to one over 
a wide range of frequencies. Typical upper limits of the power-law decaying 
range are from 100 to 10 4 Hz, and a typical lower limit is 10 _2 Hz. (In some case 
the lower limit has been reported to be extended down to 10 _7 Hz without ob- 
serving the shoulder frequency.) [22] Until now, such beautiful 1/ / fluctuations 
have never been reproduced from any Hamiltonian systems with many degrees 
of freedom. In case of one-dimensional FPU /3-model, it was found [23], for 
large N, that the shoulder frequency fo of the lowest mode energy spectrum 
decreases with increasing N, while the spectrum is close to Lorentzian: a ~ 2 
in a frequency range / > f . In other words, the energy fluctuations of the 
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FPU /3-model cannot be said as l//-type even if the time scale ~ l//o might 
be long. In contrast, the spectra for diatomic gas obtained here deviate from 
Lorentzian and have an exponent a smaller than two. 

Brillouin scattering experiments on quartz [24] showed that phonon number 
fluctuations (energy fluctuations of each normal mode) exhibit l/f a spectra 
over a frequency range from 10~ 4 to 1CT 1 Hz. In similar experiments on liquid 
water [25], 1/ f a spectra were also observed, in which the frequency range is 
even wider: 1CT 4 to 10 1 Hz. In both cases, a is close to one and no shoulder 
appears during the observation time which is of orders of hours (~ 10 4 s). As 
stated above, one-dimensional FPU models are unsatisfactory to explain such 
1/f noise. On the other hand, modified FPU models having internal degrees 
of freedom seem to be a better candidate. Actually, the analysis of a modified 
one-dimensional FPU model [6] is suggestive, in which the long time scale was 
proved to remain in the thermodynamic limit. At this stage, we think much 
of investigation into the power spectra of such systems to know whether the 
power-law decaying range could be wide enough in the thermodynamic limit 
and, more importantly, whether a could approach one owing to the internal 
degrees of freedom. 
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